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Recent fully nonlinear, kinetic three-dimensional simulations of magnetic reconnection [l|] evolve 
structures and exhibit dynamics on multiple scales, in a manner reminiscent of turbulence. These 
simulations of reconnection are among the first to be performed at sufficient spatio-temporal res- 
olution to allow formal quantitative analysis of statistical scaling which we present here. We find 
that the magnetic field fiuct nations generated by reconnection are anisotropic, have non-trivial spa- 
tial correlation and exhibit the hallmarks of finite range fluid turbulence; they have non-Gaussian 
distributions, exhibit Extended Self-Similarity in their scaling and are spatially multifractal. Fur- 
thermore, we find that the field J • E is also multifractal, so that magnetic energy is converted to 
plasma kinetic energy in a manner that is spatially intermittent. This suggests that dissipation 
in this sense in collisionless reconnection on kinetic scales has an analogue in fiuid-like turbulent 
phenomenology, in that it proceeds via multifractal structures generated by an intermittent cascade. 

PACS numbers: 94.30.cp 94.05. Lk 



Magnetic reconnection is a fundamental process that 
converts magnetic energy into various forms of plasma ki- 
netic energy. It is thought to occur in a variety of space, 
astrophysical and laboratory applications, with parame- 
ter regimes spanning from collisional to highly collision- 
less plasmas (e.g. see Ref. [H and references therein). 
While many studies have focused on laminar initial con- 
ditions, it is now widely recognized that the influence of 
turbulence remains a major uncertainty. Depending on 
the application, the turbulence may arise from a spec- 
trum of instabilities within the reconnection layer or from 
pre-existing magnetic fluctuations in the ambient plasma. 
Within the magnetohydrodynamic (MHD) model, there 
has been progress on both ideas - either by starting with 
a laminar current sheet to explore instabilities [3|-l6| or 
by directly driving turbulence (7l-[ll|. 



Moving beyond the MHD model into kinetic regimes, 
most research has focused on initially laminar current 
sheets within a variety of descriptions [12] including two 
fluid, hybrid and fully kinetic simulations, which allow 
a complete description of the electron physics respon- 
sible for breaking the frozen-flux condition in collision- 
less parameter regimes [l3|, [l^ . As larger kinetic simula- 
tions became possible, one surprising result was that the 
nonlinear evolution of reconnection produced extended 
electron-scale current sheets, with half-thickness on the 
order of the electron inertial length and leng ths that can 
extend beyond the ion inertial scale |15l4l9|. These pre- 
dictions have since been confirmed in spacecraft obser- 
vations [ioj. While the precise details depend on the 
strength of the guide field (i.e. magnetic shear angle), it 
has been demonstrated that electron pressure anisotropy 
plays a key role in setting up and driving these layers 



|2l| . Thus the existence of these structures is now well 
accepted and a variety of two-dimensional (2D) kinetic 
simulations have shown the layers can become unstable 
to secondary magnetic islands 15-18| leading to time 
dependent scenario. Recent extensions of these kinetic 
simulations to 3D have demonstrated that the tearing 
instability within these electron layers has much greater 
freedom to develop and gives rise to numerous magnetic 
flux ropes [1]. The subsequent nonlinear interaction of 
these flux ropes is seen in the simulations to lead to the 
self-generation of structures on multiple scales within the 
initially laminar ion-scale current layer. The question 
then arises - is the multi-scale nature of these flux ropes 
important for dissipation in the sense of conversion of 
magnetic to kinetic energy? The key properties of turbu- 
lence in this context are that it cascades energy in a scale 
invariant and intermittent manner. Thus turbulence pro- 
vides a mechanism to form a spatial multifractal fleld of 
coherent structures across a broad range of spatial scales 
which can then contribute to dissipation of magnetic fleld 
energy into plasma kinetic energy. 

In this Letter we show that these fully kinetic 3D PIC 
simulations of magnetic reconnection [l| do indeed ex- 
hibit the hallmarks of intermittent turbulence. A key 
property of turbulence is that it can be characterized and 
quantifled in a robust and reproducible way in terms of 
ensemble averaged statistical properties of fluctuations. 
In an inflnite medium, fully developed turbulence ex- 
hibits statistical scale invariance in fluctuations in the 
bulk quantities that describe the flow. Either when 
the turbulence is not completely evolved (low Reynolds 
number), or the system is of flnite size, a generalized 
scale invariance or extended self similarity (ESS) 0] still 
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holds in both hydro dynamic and magnet ohydro dynamic 
(MHD) turbulent flows, as seen for example, in the solar 
wind [e.g., [23l - [25| and in a solar quiescent prominence 
(26| . A key discriminator of turbulence is that fluctu- 
ations have non Gaussian probability density functions 
(PDFs) with moments that exhibit multifractal scaling. 
Remarkably, we find that the magnetic field fluctuations 
on kinetic scales in these 3D simulations of reconnection 
have this quantitative statistical character, specifically, 
they show non-Gaussian statistics and ESS consistent 
with multifractal scaling. Furthermore, dissipation, in 
the sense of conversion of magnetic to kinetic energy, oc- 
curs on a multifractal field. This suggests that dissipa- 
tion in collisionless reconnection on kinetic scales has an 
analogue in fluid-like turbulent phenomenology, in that 
it proceeds via multifractal structures generated by an 
intermittent cascade. 

We focus on 3D simulations of undriven magnetic re- 
connection in collisionless plasmas [1]. These petascale 
simulations use the kinetic particle-in-cell code VPIC, 
which solves the full set of relativistic Maxwell's equa- 
tions. The simulation is initialized with a Harris current 
sheet, with magnetic field B = Bxot(inh{z / \)ey^^ Byo^y^ 
where ex and Gy are unit vectors, Byo = B^o is a uni- 
form guide field and \ = di \s the initial half- thickness, 
where di is an ion inertial length. Here the ion elec- 
tron mass ratio is mi/rrie = 100, which implies that 
di is lOde, the electron inertial length. The simula- 
tion employs open boundary conditions in the x and 
z directions and periodic boundary conditions in the y 
direction. The dimensions of the simulation grid are 
LxXLyXLz = 2, 048 x 2, 048 x 1, 024 cells, corresponding 
to 70di X 70di X 35di. We consider a slice of the 3D simu- 
lation grid in the X-Z plane at Y = 3bdi and at the time 
tQci = 78 of the simulation, where Qci = eBxo/{mic). 
This time slice corresponds to a middle phase of the mag- 
netic reconnection in which the turbulence quantitatively 
reaches its peak of evolution as seen in the power level of 
the region of scaling in the power spectral density (PSD). 
The overall time evolution of the simulation is shown in 
[H . We also analysed an early and late phase of the pro- 
cess at t^ci = 40 and 98 respectively, and find that the 
reconnection generated fluctuations evolve in time: at 
early time the power in these fluctuations is weak though 
above the noise, then it grows in amplitude up to the mid- 
dle phase where the total integrated field energy density 
over the turbulent region of the power spectral density is 
~ 0.1 that of the background field. The reconnection rate 
in these kinetic simulations is fast, Vin/VA ^ 0.1, as is 
the case for kinetic simulations of reconnection where the 
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flow is almost laminar [e.g. |27|, and references therein]; 
hence these kinetic simulations do not show correlation 
between the level of turbulence fluctuations and the re- 
connection rate as has been postulated for reconnection 
on MHD scales ji,[28|. 

Figure [T] shows the 2;-component of the magnetic field. 



FIG. 1: Bz in the X-Z plane at Y=35di and at the time 
tQci — 78 of the simulation. Solid black lines are the cuts 
chosen for the analysis of the magnetic fluctuations, while 
squares indicate the regions over which we perform the box 
counting analysis of J • E. 



Bz^ at tflci = 78 corresponding to the time-space slice 
for which we present analysis here. We take cuts on the 
simulation grid along the horizontal direction and label 
them as cut 1 to 6 (see Fig. [1]). Cuts 2-5 are within 
the reconnection region where the field is strongly fluc- 
tuating, while cut 1 and 6 lie where there is negligible 
signature of reconnection. The latter are taken to quan- 
tify the effects of the PIC simulation noise. Here we 
focus upon the analysis of the z-component of the mag- 
netic field. The reconnection generated fluctuations are 
highly anisotropic in character and have clearest signa- 
ture in their z-component which is perpendicular to the 
x-y plane of the macroscopic field of these simulations. 
This anisotropy parallels what has been recently observed 
both in kinetic range turbulence in the solar wind [e.g., 
2^, 30| and in a reconnection jet jsif. 

Fluctuations associated with fully developed intermit- 
tent turbulent flows are characterized by non- Gaussian 
probability distribution functions (PDFs). Under the 
assumption of statistical stationarity and homogeneity, 
fluctuations of a field / on length scale L along a given 
direction r are defined as SI{L) = I{r-\-L) — I{r). Figure 
[2] shows the PDFs of the fluctuations SBz{L) at different 
scales L for cut 4 (left panel) and cut 1 (right panel). We 
can see that the fluctuations of cut 4 (reconnection gener- 
ated fluctuations) follow a non- Gaussian distribution at 
scales L in the range Ade < L < 25(ie, while fluctuations 
associated with cut 1 (noise) are Gaussian distributed 
up to scales L ~ lOOdg = lOd^. The latter suggest that 
the PIC noise behaves as Brownian noise. We recover 
consistent statistics also for the other cuts considered, 
that is, Gaussian distributions for cut 6 at scales L up to 
~ lOdi and non- Gaussian PDFs for cuts 2,3 and 5 within 
the range Ade < L < 25de. We identify the latter as a 
potential range of turbulence, which we now test. 

A central characteristic of turbulence is that the PDFs 
of non- Gaussian fluctuations at different scales are re- 
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FIG. 2: PDFs of the fluctuations ^B^L) = B,{X+L)-B,{X) 
of cut 4 (left panel) and cut 1 (right panel) for six different 
values of the space lag L. In both panels, all PDFs are cen- 
tred on the mean value < 6Bz{L) > and normalized to the 
standard deviation a to allow comparisons with a Gaussian 
distribution (dashed black lines). All curves are shifted along 
the ^/-direction for clarity. 
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FIG. 3: Log- log plots of 5*2 against S3 for all the cuts of in- 
terest. Dashed black lines correspond to the linear regression 
fits in the potentially turbulence range Me < L < 25de for 
cuts 2 to 5 and within the range Idi < L < lOdi for cuts 1 
and 6. All curves are shifted along the ^-direction for clarity. 



lated by a multifractal similarity or scaling. We now test 
for this by examining the generalized structure functions 
(GSFs) of the magnetic field fluctuations SBz{L)^ defined 
as Sp{L) =< \SBz{L)\P >, where the angular brackets in- 
dicate an ensemble average over r, implying approximate 
statistical homogeneity. In infinite range, fully devel- 
oped turbulence, one would expect the GSFs to scale as 
Sp{L) ~ L^^P\ where the ({p) are the scaling exponents. 
For turbulence in a finite domain, or for turbulence that 
is not fully developed, a generalized scale invariance or 
ESS can hold, with Sp{L) - G{L)^^p\ where the func- 
tion G{L) is an initially unknown function that depends 
on the largest scale physical structures [H, [H, [111 . 
While for fractal fields ({p) is linear in p, intermittent 
turbulence is realized by multifractal topology (e.g. 
and the ({p) are non linear in p. In finite range turbu- 
lence, one does not have direct access to the scaling ex- 
ponents C{p)^ however it possible to obtain their ratios, 
C{p)/C{q)^ by plotting one structure function of order p 
against another structure function of order q. Thus, ESS 
establishes the following scaling for the structure func- 
tions 5p(L) = [^,(L)]^^^^/^^^^ [2^]. 

In Figure [3] we plot ^2 versus ^3 on a logarithmic scale 
for all the cuts considered. We see that ESS holds for 
both cuts 1 and 6, which are simulation noise, and cuts 
2 to 5, which are within the potentially turbulent range 
Ade < L < 2bde. This implies that both noise and re- 
connection generated fluctuations have a range of scale 
invariance. Now, the question that immediately arises 
as to whether these fluctuations are multifractal or not. 
We calculate all possible combinations of the scaling ex- 
ponents ratio C{p)K{q) P->Q. = 1^2,3 and 4 by plot- 
ting log(5'p) versus log(5'g) and by reading the gradients 
of the linear fits to these curves within the potentially 
turbulent range, Ade < L < 25de^ for cuts 2 to 5 and 
in the noise range, Idi < L < 10(ii, for cuts 1 and 6. 



The panels in Figure |4] show the ratios C{p)/C{q) v^^- 
sus p for g = 1 up to 4 for cuts 1 and 6 (noise) and 
cuts 2 to 5 (reconnection- fluctuations). The noise cuts 
(blue rectangles) show a linear behaviour of ({p) with p, 
consistent with a fractal field. The PIC noise thus gen- 
erates a spatial field of magnetic fluctuations which is a 
self-affine Brownian noise, which shows fractal scaling. 
Importantly, it is clearly distinguishable from the recon- 
nection generated fluctuations of cuts 2 to 5 (green dots), 
which instead consistently show a non linear dependence 
of ({p) on p within uncertainty. The reconnection gen- 
erated structures are thus robustly characterized by a 
multifractal spatial field. At minimum, this suggests a 
new signature of reconnection outflow regions that could 
provide a method for observational identification, as has 
indeed been recently observed [lH. However, this is also 
a key signature of a multifractal intermittent turbulence 
phenomenology. It suggests that dissipation, in the sense 
of conversion of magnetic to kinetic energy, in collision- 
less reconnection on kinetic scales has an analogue in dis- 
sipation in fluid-like turbulent phenomenology, in that it 
proceeds via a spatial multifractal field of structures gen- 
erated by an intermittent cascade. If this is the case, then 
the spatial dissipation field will also be multifractal. 

We now test this idea by directly quantifying the spa- 
tial topology of the field J • E. We show results for the 
three components of J • E, at the same time-space sHce 
of the simulation discussed above. We perform the clas- 
sical box counting method to calculate both global and 
local fractal dimension [35] for a reconnection dominated 
turbulent region and a region where noise is dominant. 
The box-counting method consists of dividing a spatial 
region into boxes of size L using a regular grid, and then 
counting the number of boxes N{L) that contain non- 
zero values of a discretized spatial field. We consider the 
magnitude of each component of J • E and the discretized 



4 





Moment p 



Moment p 



Box size, L [d ] 



FIG. 4: C{p)/C{q) versus p for cuts 1 and 6 (blue rectangles) 
and cuts 2 to 5 (green dots) and for ^ = 1 (top- left), 2 (top- 
right), 3 (bottom- left) and 4 (bottom-right). Solid and dashed 
red lines are the best fits to the curves within uncertainty for 
the turbulent (cuts 2-5) and noise (cuts 1,6) cuts respectively. 



FIG. 5: Box-counting method. Plot of the number of 
boxes, N(L), versus the box size L within a reconnection dom- 
inated turbulent region (red square in Fig[T]) and a noise re- 
gion (yellow square in Fig|T]) for each component of J -E. The 
inset figure shows the corresponding local dimension, n(L) = 
d ln(N(L))/d ln(L), against the box size L. 



values are non-zero where it exceeds a threshold, we test 
the robustness of our results by varying this threshold. 
For sufficiently small thresholds, this yields the topology 
of the noise field, but for thresholds above the noise, we 
obtain the topology of the turbulence. Box counting the 
thresholded turbulent field then gives its spatial topology 
in the absence of noise, without the need of a filtering or 
averaging. This method can thus probe spatial struc- 
tures on length scales where the noise power dominates 
the overall signal power provided that at least in some 
locations, the signal is above the noise threshold. For 
fractal geometries N(L) is expected to depend linearly 
upon L, while non-linear trends of N(L) against L indi- 
cate a multifractal field. In Figure [5] we plot N(L) versus 
L within the reconnection dominated turbulent region 
and the noise region indicated respectively by the red 
and yellow squares in Figure [TJ The noise region shows a 
linear behaviour of N{L) with L indicating that the PIC 
noise is fractal, on the contrary, in the turbulent region 
the plot introduces curvature for scales L smaller than 
^ 25de. The inset of Figure [5] also shows how the local 
dimension, n(L) = d ln(N(L)) / d ln(L), varies with the 
scale L. Within the noise region, the fractal dimension 
is roughly constant as L varies, while it changes with 
the scale L in the turbulent region again, a signature of 
multifractality. Thus dissipation, in the sense of energy 
transfer to the plasma via J • E, occurs in a spatially 
intermittent manner. 

In conclusion, recent fully kinetic (PIC) simulations in 
3D reveal that reconnection is dominated by magnetic 
structures on multiple scales which manifest highly vari- 
able dynamics suggestive of turbulence. We have quan- 
tified the ensemble averaged statistical properties of the 



spatial fields of fluctuations in the magnetic field and in 
energy transfer to the plasma. The magnetic field fluctu- 
ations are anisotropic and exhibit the hallmarks of finite 
range fluid turbulence; they have non- Gaussian distri- 
butions, exhibit Extended Self- Similarity in their scal- 
ing and are spatially multifractal. These signatures are 
recovered quite robustly across the regions in the simu- 
lation domain where reconnection is actively generating 
fluctuations. This potentially offers a new observational 
test for reconnection regions using in-situ observations, 
so that for example recent observations of non- Gaussian 
fluctuations in a turbulent jet [31] (see also [sgI, [37| ) 
could be tested for ESS and non-linear ratios of expo- 
nents as found here. Furthermore, the fact that we also 
find that the spatial field of J • E is multifractal sug- 
gests that the turbulence converts some of the magnetic 
energy to plasma kinetic energy in a spatially intermit- 
tent manner. This suggests that dissipation in the sense 
of energy transfer to the plasma in collisionless recon- 
nection on kinetic scales has an analogue in fluid-like 
turbulent phenomenology, in that it proceeds via multi- 
fractal structures generated by an intermittent cascade. 
This provides a starting point for theoretical models of 
heating in collisionless reconnection; it also suggests that 
existing analytical and quasi-analytical models of recon- 
nection that do not take into account the development of 
turbulence in the reconnection layer may need revision. 
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